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Abstract. This paper is devoted to the numerical approximation of a nonlinear tempera- 
ture balance equation, which describes the heat evolution of a magnetically confined plasma 
in the edge region of a tokamak. The nonlinearity implies some numerical difficulties, in 
particular long time behavior, when solved with standard methods. An efficient numerical 
scheme is presented in this paper, based on a combination of a directional splitting scheme 
and the IMEX scheme introduced in [1]. 
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1. Introduction 

The description and simulation of the transport, especially the turbulence of magnetically 
confined fusion plasmas in the edge region called scrape off layer (SOL) of a tokamak, is 
nowadays one of the main problems for fusion generated energy production (ITER). The 
understanding of the physics in this edge region is fundamental for the performances of the 
tokamak, in particular the plasma-wall interactions as well as the occurring turbulence have 
an important impact on the confinement properties of the plasma. From a numerical point 
of view, an accurate approximation of the plasma evolution in the edge region is essential 
since energy fluxes as well as particle fluxes at the boundary are used as boundary conditions 
for the mathematical model applied to describe the plasma evolution in the center region 
(core) of the tokamak. The physical properties of these two regions (core/edge) are rather 
different, so that different models are used for the respective plasma-evolution modeling: the 
gyrokinetic approach for the collisionless core-plasma and the fluid approach for the colli- 
sional edge-plasma. 

A large variety of models can be found in literature [5, 9] for the description of the SOL, 
based on various assumptions and aimed to describe different physical phenomena. We shall 
concentrate in this paper on the T0KAM3D model, introduced in [s]. The aim of this model 
is the investigation of the instabilities occurring in this plasma edge region, as for example 
the Kevin-Helmholtz instability, the electron-temperature-gradient (ETG), ion-temperature- 
gradient instabilities (ITG) , etc. 

The T0KAM3D model is based on a two-fluid description (ions, electrons) and consists of 
the usual continuity equation, equation of motion and energy balance equation, closed by the 
so-called "Braginskii closure". These equations are 

dtUa + V • {riaUa) = Sna , 



(1.1) 



^jHa [dtTa + {Ua ■ S/)Ta] + PaV • n« = -V ■ Qa - ■ Vu^ + Qc 



2 "'a I'-'t'^a I \ '^a v jj^ai i Fa ^ "•a — ^ Ha '-'-a ■ ^ ^a « 5 

where Ua is the particle density (a = e for electrons and a = i for ions), Ua the velocity, 
Ta '■= UaUa the particle flux, the particle mass, the particle charge (eg = —1 for 
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electrons and = 1 for ions), pa the pressure, lia the stress (viscosity) tensor, Sna a 
particle source term (coming from the core plasma), Ra the friction force due to collisions, 
Tq the temperature, the energy flux and finally Qq, the particle exchange energy term, 
due to collisions. In the Braginskii closure, the pressure is specified as pa '■= naTa (perfect 
gas assumption), the plasma viscosity is supposed negligible, such that V • Ha = and 
n : Vuq, = and the energy flux qa is supposed to have a diffusive form, given in terms of 
the temperature gradient, as follows qa ■= —Ka'VTa (coming from the Fourier law) with 
the thermal conductivity coefficient. The energy exchange term Qa is taken under the form 



Qa := ±3^^(re-r,), 



where Te is the electron-ion collision time. 

Due to the high complexity of the problem, several other hypothesis are assumed, permitting 
to concentrate on the desired features and to filter out the insignificant/disturbing details. 
These hypothesis, as for example the quasi-neutrality rig ~ rij, are not detailed here and we 
refer the reader to the more physical works [3, 9, 10]. 

Several difficulties arise when trying to solve numerically the system (1.1). We shall con- 
centrate in this paper only on the temperature equation, which requires at the moment still a 
lot of effort, due to its inherent numerical burden. The resolution of the two other equations 
was the aim of the PhD thesis [8]. The numerical difficulties in solving the temperature 
equation are firstly related to the thermal conductivity coefficients, which depend on the 
temperature itself, leading thus to a non-linear problem. Secondly, the strong magnetic field 
which confines the tokamak plasma introduces a sharp anisotropy into the problem. Indeed, 
the charged particles gyrate around the magnetic field lines, moving thus freely along the 
field lines, but their dynamics in the perpendicular directions is rather restricted. Quantities 
as for example the resistivity or the conductivity, differ thus in several orders of magnitude 
when regarded in the parallel or perpendicular directions. Finally, boundary conditions have 
to be imposed, which is a rather delicate task from a physical, mathematical and numerical 
point of view. 

Let us now present in more details the model we are interested in. In this paper, we shall 
study a simplified version of the temperature evolution equation, which contains however all 
the numerical difficulties of this last one. We shall focus on how to handle with the nonlinear 
terms and the boundary conditions, the high anisotropy being the aim of a forthcoming work 
[1, 6]. The simulation domain = (0,1) x (0,1) with boundary is presented in Figure 
1. It consists of a periodic core region, separated by a Separatrix from the non-periodic 
SOL region. Its axes represent the direction parallel to the magnetic field lines (s) and the 
radial direction (r). We assume in this paper that all quantities are invariant with respect 

5/2 

to the poloidal angle ip. The parallel thermal conductivities Kg || depend on Ta whereas the 
perpendicular ones Ks,±, governed by the turbulence, are independent of the temperature [2]. 

The system we are interested in, is composed of the evolution equation 



(1.2) dtTa-ds{Ki^,aT^/^dsTa)-driK^,adrTa) = ±f3a{Te-Ti), (s, r) G , 
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completed with the boundary conditions 

r drT^ = -Q^,„ , r = , s G (0, 1) , 



(1.3) 



drTa = ^, r = l, sG (0,1), 
K\\^^T^'^ dsT^ = 7aTo„ r G (1/2, 1) , s = , 
i^,i,„Tj/2 5,T« = -7„r„, re (1/2,1), s = l 



t T{t,0,r)=T{t,l,r), rG (0,1/2), 
and the initial condition 

(1.4) TaiO) = T^,o > . 

The diffusion parameters < K±^a ^ ^\\,a core- heat flux > are considered 

as given. The non-linear boundary conditions at the limiter express the fact, that we have 
continuity of the heat fluxes at the boundary. Indeed, the heat flux q := 7r||T at the 
boundary is given as the sum of a diffusive and a convective term, like 

5 

-fr\\T = -T\\T - \K^^\dsT , r||=n-U||. 

At s = the particle velocity < is negative, whereas at s = 1 we have > 0, which 
gives rise to the boundary conditions in (1.3). The constant ja is different for electrons and 
ions, in particular 7i ~ for ions and 7e ~ 5/2 for electrons. In the case of ions, we have 
thus homogeneous Neumann boundary conditions at the limiter. 
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Figure 1 . The 2D domain. 



The outline of this paper is the following. In Section 2, we will focus on the ID nonlinear 
parabolic problem 

dtT-ds{KilT^^^dsT) = 0, 

completed with the nonlinear boundary conditions in s = 0, 1. A mathematical study is firstly 
performed. Then, explicit, implicit and IMEX schemes are compared for the resolution of 
this ID problem, with respect to precision and simulation time. In Section 3 we consider the 
complete 2D problem for one species (without the source term). A directional Lie splitting 
method is used in order to transform the 2D problem in two ID problems and to apply the 
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results of the previous section. Finally, in Section 4 we solve the complete 2D ion-electron 
coupled problem. The shapes of the different electron/ion temperatures are compared. 

2. The ID nonlinear problem 

Let us consider in this section the ID nonlinear problem, corresponding to the temperature 
balance equation in the parallel direction, i.e. 

dtT - ds{K. !r|5/2a,r) = o, (t, s) g m+ x (o, i), 



(2.1) 



iC|i|r|^/2^,r = 7T, s = o, 

KATf'^dsT = --iT, s = 1, 



r(o,-) = to, 

where 7 > is a given constant, K\\ G L~(f]), r° G ^^(f]), with r° > and i^Tn > almost 
everywhere. Let us denote in this section the domain by J7 = (0, 1) and the time-space 
cylinder by Q := M"*" x (0, 1). The aim of this section is to study from a mathematical point 
of view this equation and to introduce an efficient numerical scheme for its resolution. From 
a physical point of view, problem (2.1) describes the rapid diffusion process of the initial 
temperature and the outflow through the boundary. 

2.1. Mathematical study. Before starting with the numerical discretization, we first es- 
tablish some properties of the ID diffusion problem (2.1), like existence, uniqueness of a 
solution, positivity etc. To simplify the presentation, we shall assume for the present study 
that K\\ = 1, the general case being treated equally. 

We also denote p > 2 and p' its conjugate number 1 < p' := < 2. The diffusion 
coefficient can now be written as a(T) := |rp~'^. Moreover, let us define the primitive 

A(r) := / a(x)dx = \T\P-^T. 

Jo P - 1 

With these notations, the diffusion equation can be simply rewritten under one of the two 
forms 

dtT - ds {a{T)dsT) = , dtT- dss (A(T)) = . 

We shall now introduce the concept of weak solution of problem (2.1) and state the exis- 
tence/uniqueness theorem. 

Definition 2.1. Let us consider T^ G L^(i7) and define W C L^iQ) as the space 



:= {t G L\Q), ITI'-^T G LL(M+, i/i(17)), dtT G LfjR+, (t^^'^ (f)))*)} 



and we denote by T> = C^(IR+, VF"'^'^(il)) the space of test functions. Then the temperature 
T &yV is a weak solution to (2.1) if and only it satisfies 

[ [ T{t,s)dtip{t,s)dsdt- [ [ \T{t,s)\P-^dsT{t,s)dsip{t,s)dsdt 
Jr+ Jn Jr+ Jn 

-7/ [T{t,l)^{t,l) +T{t,0)ip{t,0)]dt + [ T^{s)(f{0,s)ds = 0, G P. 
Jr+ Jn 

Remark that all the terms in this variational formulation are well-defined. Moreover, we 

observe that a function T G LP(M+ x Q) satisfying dtT G Lf^^(M+, (T^^'f (0))*) belongs to 

C{[0,t], {W^'P{^1))*), for all r > by Aubin's Lemma, such that the initial condition is 

well-defined. 



NUMERICAL STUDY OF A NONLINEAR HEAT EQUATION FOR PLASMAS 5 

Theorem 2.2. Let £ L^(il) with > 0. Then, there exists a unique weak solution 
T €W of (2.1), which satisfies T e L°°(M+, L2(J^)), T > almost everywhere and 

|lim-)iii.(n)<o. 

The proof of this theorem is decomposed in several steps. For the begimiing, we shall 
suppose that € L^{Q), with HT'^Hoo < M and fixed M > 0. A truncation can be done, 
for more general G L^{U). 

Two main difficulties arise in the mathematical study of (2.1), the nonlinearity and the 
degeneracy, which means that the equation changes its type there where T = 0. 

Proof. We shall first regularize the problem, in order to avoid the degeneracy. Then, in a 
second step, we shall treat the nonlinearity via a fixed point argument. Finally, a priori 
estimates shall help us to pass to the limit, in order to deal with the degenerate problem. 
Let us thus detail these steps. 

First step: Regularization. Let < e < 1 be fixed and let us define the regularized 
diffusion coefficient 

a,,Af(T) := [e2 + min(|T|2,M2)]' 
and the corresponding primitive 

rT 



p-2 
2 



K,m{t) ■■= [ 

Jo 



0'e,M{x) dx . 



The diffusion coefficients being now bounded from below and above, standard arguments 
allow to prove that the regularized problem 

r dtT.^M - ds{a,,M{Te,M)dsT,,M) = 0, (t, s) G M+ X (0, 1), 



(2.2) 



ae,M{Te,M)OsT^,M = 7^e,M, S = 0, 
<le,M{Te,M)dsTe,M = —jT^M' ^ = 1, 



, r(o,-) = ro, 

has a unique weak solution T^^m £ L'^ (M.'^ , (Q,)) such that it satisfies the following varia- 
tional formulation: for any ip £ C^iW'^ , H^{Q,)) 



(2.3) / / T^^Mit,s)dtip{t,s)dsdt - / ae,MiTe,M)dsT^,Mit, s)dsip{t, s) ds dt 
-7 / [T,^Mit,l)HtA)+Te,Mit,0)Ht,0)]dt + [ r°(s)99(0,s)ds = 0. 



These arguments are based on the Schauder fixed point theorem, applied on the mapping T 
: Bpi I— 7> Bpt with 

Br := {v G L\Q), \\v\\l2(^q) < R} , 

where for v G Br we associate Tv the solution of the linearized problem associated to (2.2). 
2nd step: a priori estimates. In order to pass to the limit e — )• 0, we will need some a 
priori estimates for the solution T^^m, independent of e. Taking in the variational formulation 
(2.3) as test function T^^m, yields first 

J [ \T,,M{t,s)fds + [ [ a,^M{Te,M)\dsTe,M\^dsdT 

^ Jn Jo Jn 



+ 7 / [|r,,M(T,l)P + |r,,Af(r,0)p]dr = \ f \T' 
Jo ^ Jn 



Ol „\\2 



ds 
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which imphes that for all t > 

ll^e,Af(*)||L2(f^) < ||T°||i2(n), 



/ / a^^M{T^^M)\dsT^,M\'^dsdT < ||T°|||2(m • 
K Jo Jn 



p — 2 

This shows also, that the sequence {\T^,M\~T,^M}e is bounded in L^{R+,H^{n)) and hence 
{Te,M}e bounded in U'(Q). Moreover by standard arguments for parabolic problems we 
deduce than, that {dtT.^M}, is bounded in Lp'(R+, {W^'P{n))*). 

Third step: passing to the limit. The a priori estimates of the last step permit us to 
show, that there is a sub-sequence and a function Tm G L'^{Q), such that 

r,,M ^ Tm in L^{R+ x Q) as e ^ 0. 

Moreover, from standard compactness arguments [7], we show that the following set of mea- 
surable functions 

T := \u:m+ ^ {w^'P{n)y, [ [ \u\p-^\dsu\^dsdt<c , dtueLP'{R-^,{w^'P{n)y)] , 

[ Jr+ Jn J 

is compactly embedded in x $7), implying thus that, up to a sub-sequence 

Te,M Tm, in LP(M+ x n) , as e ^ 

p—2 

and then T^^m Tm a.e. in Q when e goes to zero. Furthermore, since {\T^,M\^^Tf:^M}e is 
bounded in L2(M+, H^{n)), one has 

{TeMl'^T.^M ^ {TmI'^Tm , in L\R+,H\n)), as e ^ 0, 



implying by the weak continuity of the trace application 

\T,^m\^T,^m \Tm\^Tm, in l2(R+ x dn), as e ^ 0. 
Finally, we also have using the same arguments, when e ^ 

KM{Te,M) ^ ^m{Tm), in L\R+,H\n)), 

dtT.^M ^ dtTM, in LP'(M+, (T^^'^ (J]))*). 

All these convergences permit us now to pass to the limit in the variational formulation (2.3) 
in order to show the existence of a weak solution of problem (2.1). This solution is even 
unique and satisfies the maximum principle, which can be shown as in step 2. □ 

2.2. A finite volume approximation. In this section, we propose to derive a numerical 
scheme for (2.1) in which we apply a finite volume approach for the discretization in the space 
variable. Let us consider a set of points (si-i/2)o<i<ns of the interval (0, 1) with s_i/2 = 0, 
*ns-i/2 = 1 S'Hd rig + 1 represents the number of discrete points. For < z < — 1, we 
define the control cell Ci by the space interval Cj = (sj_i/2) •Si+i/2)- We also denote by Si the 
middle of Ci and by Asj the space step Asj = Sj+1/2 — ■Si-1/2 where we suppose that there 
exists ^ € (0, 1) such that 

(2.4) ^As < Asi < As, Vi G {0, ...,ns - 1}, 

with As = maxj As,. 

We shall construct a set of approximations Ti{t) of the average of the solution to (2.1) on 
the control volume Ci and first set 

T^ = ^ I Us)ds. 
^Si Jc^ 
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Applying a finite volume discretization to (2.f), Tj is solution to a system of ODEs, which 
can be written as 



(2.5) 



-7r(*) = X ' 0<i<ns-l, 

at Asj 



Ti{t = 0)=T^, < i < n, - 1, 
where the numerical flux is given by 

2.6 -^^+1/2 = ^ ^ ^^ , i = 0,...,n,-2. 

7 Asj+i + Asi 

Moreover, at the boundary s = and s = 1, we apply the boundary conditions, 

{+7 To, ifz = -l, 
-■jTn^-i, ifi = ns-l. 

Note that the above discretization on space is first order due to the loss of precision at the 
boundary. To complete the discretization to the system (2.f ), the finite volume scheme (2.5)- 
(2.7) has to be supplemented with a stable and consistent time discretization step. In the 
following we present different time discretizations starting from classical explicit and implicit 
schemes and then propose a stable and accurate numerical approximation. 

2.3. Time explicit discretization. We denote by At > the time step, = nAt for any 
n S N and T" is an approximation of the solution T to (2.1) at time t". Then, we apply a 
backward Euler scheme to (2.5)-(2.7), which yields 



(2.8) 



At Asi 
Tf = ToA, < i < n, - 1, 



with -^7+1/2 ^^-^ (2-6)-(2.7) computed from the approximation at time T". 

Classically, to guarantee the stability of the scheme (2.8), the time step At is restricted by 
a CFL condition. 

Proposition 2.3. Consider that the initial datum Tq is nonnegative and Tq G L°°(0, 1) and 
assume the stability condition 

(2.9) At < ^ 



max ( ^^||ro||^^,7As 



where ^ is given in (2.4)- Then, the numerical solution {T^)i^n obtained by the explicit scheme 
(2.8) is stable and converges to the exact solution to (2.1). 

We don't give the proof of this result since it is similar to the proof of Proposition 2.5 
presented in the next section. Unfortunately, this simple scheme is not really efficient since 
it becomes costly when the mesh is very fine, the constraint on the time step becoming too 
restrictive. 



2.4. Time implicit discretization. To avoid the restrictive constraint on the time step 
(2.9), an implicit scheme is more suitable. Therefore, we consider the finite volume scheme 
(2.5)-(2.7) to the system of equations (2.1), but apply a forward Euler time discretization. 
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This yields, 
(2.10) 



i-1/2 



At 

TP =Tr 



0,41 



Asi 

< i < — 1, 



< i < n. 



1, 



with .7^^1/2 ^^-^ (2-6) computed from the approximation at time T"""*"^. Hence, a fully 
nonlinear system has to be solved at each time step. 

The scheme (2.10) coupled with (2.6)-(2.7) is uniformly stable and leads to a numerical 
approximation which converges to the exact solution to (2.1). 

Theorem 2.4. Consider that the initial datum Tq is nonnegative and Tq € L°°(0, 1). Then 
the numerical solution given by the implicit scheme (2.10) coupled with (2.6)- (2.7) is uni- 
formly stable in L°°(M"^ x (0,1)) and converges to the weak solution T of (2.1) when h = 
(At, As) goes to zero. 

We start with a stability result and then prove convergence of the numerical solution to 
the unique weak solution by consistency of the scheme. 

Let us first investigate the stability property and prove some a priori estimates on the 
numerical solution uniformly with respect to the mesh size h. 

Proposition 2.5. Consider that the initial datum Tq is nonnegative and Tq € L°°(0,1). 
Then the numerical solution given by the implicit scheme (2.10) coupled with (2.6)-(2.7) is 
unconditionally stable, i.e. 



(2.11) 

and 

(2.12) 



71. — 1 



n. — 1 



i=0 i=0 

Moreover, the following discrete semi-norm is uniformly bounded 



(2.13) 



Nt ns-2 
n=0 j=0 



{T! 



-n7/2 



T! 



< c, 



Asi + Asj+i 

where the constant C > only depends on the initial datum Tq. 
Proof. Let us consider a convex function G Ci(M,M), then we have 
(2.14) 0(7;"+^) - (PiTD < (t)'{Tl'+^){T^+^ - TP). 

Thus, we multiply the scheme (2.10) by At Asi (^'(T""''^) and sum over i G {0, . . . , n<j — 1}, it 
gives 



n. — 1 



ris — l 



n»-l 



^ As../<(rf+i) - Y.^'^<t>{Tr 



1=0 



i=0 



~n+l 
-1/2 



< AtY,^'{T-^')[r;^^\)^-T-_ 

1=0 

< -At J-f^Y/2 {^'(T-,f ) - <t>'{Tr')) 



i=0 



AtF-_tj,4>'{T-^^) + At j:+\/2 0'(r„"+\). 
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Using the definition of the numerical flux (2.6) and the discrete boundary conditions (2.7), 
we get 

ris — l Us — l 

^ As,0(Tr+i) - As,<P{Tn < -^AtcP'{T^+')T^+' - 7 At,/<'(r„"+\)T„"+\ 

i=0 i=0 

Observing that a similar inequality holds true when (j){x) is only Lipschitzian, we take (f){x) = 
x~ , and prove the nonnegativity of the approximation T", that is, 

< ^ A.,(Tf+i)- < As^T^r = 0. 

1=0 i=0 

Therefore, assuming that > 0, for all < « < - 1, we obtain that > for all 
< i < — 1 and n G N. Moreover, taking (/>(x) = {x — M)"*", with M = |[T''||ioo, we have 

Us-l Us — l Us — I 

< ^ Asi (1;"+^ - M)+ < - M)'^ ^ Yl " ^^)^ = 0. 

i=0 i=0 i=0 

Hence we deduce that < < M, for ah < i < - 1. 
Then we take 4>{x) = x'^/2, which yields that 

E ^ (^"^')' - E ^ (^")' <CAtY [TJX\' - T^'] ^ '^X+A. 

1=0 1=0 1 = 

and use the fact that T" is uniformly bounded to observe that 

1^7/2 _ ^7/2 1 ^ ^1^^^^ _ 

Thus, we have the following inequality 

i J;a.,(i;"+1)2 < ^E^^^(^")' 

i=0 i=0 



CAt ^ ii ^ ^ ^ ^ 



Asi + Asj+i 

Finally we sum over n € {0, . . . ,Nt} and immediately deduce that there exists a constant 
C > only depending on the initial datum Tq such that 



„ . -, Asi + Asj+i 

n=0 4=1 



< c. 

n 



2.5. Proof of Theorem 2.4. To prove the convergence of the discrete solution iT^^i^^ 
towards the weak solution T to (2.1), we construct a piecewise approximation T/^, where 
h = {At, As), such that 

Us—l 

Th{t,s) := E E Tric.is)^n^tn+i[{t), 

nSN i=0 
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From the uniform bounds proved in Proposition 2.5, we get that there exits a sub-sequence, 
stih denoted by {T^jh-, such that T/j converges to T € L°°(R+ x (0, 1)) as m ^ oo in the weak- 
* topology, whereas using (2.13) we also get that tJ^^'^ converges strongly in L^(]R+ x (0, 1) 
to Tl'\ 

Now let us prove that converges to the weak solution to (2.1) when h goes to zero. We 
consider (p G C^(M+ x (0, 1)), and we denote = ip(t^,Si). Then we multiply the scheme 
(2.10) by Asif^, and sum over z G {0, . . . , — 1} and n G N, we obtain 

4 + 4 = 0, 

with S}^ is related to the time discretization and is given by 

ris — l 

8l := a.,,(t;"+i -T;")^r, 

ngN j=0 

whereas E"^ is related to the space discretization and reads 

n.-l 



4--= AtEE 



1/2 i-1/2 

neN i=0 



'Pi 



On the one hand, we consider and perform a discrete integration by part with respect to 
rz G N. Using that if is compactly supported for large t G M+, it yields 

Us — I ris-l 



neN* i=0 i=0 



= - / I Th{t + At,s)dMt,s)dsdt-[ Th{0,s)'f{0,s)ds + < fil,ip> 
Jr+ Jo Jo 

where the additional term < fi\,ip > is given by 

<^il^> = -Y.Y.I_J / T,^dfMt,rj)dr]dsdt 



neN* i=0 
n.-l 



E / rT^ds^iO,7])drids 
i=o -^^^ 

and satisfies the following estimate 

\<^il,ip> \ < CAs [WdlMW + I|5sV5(0,.)IIli) • 
Therefore, when h tends to zero, we have 



£l^-f [ T{t,s)dt^{t,s)dsdt- [ To{s) ip{0, s)ds 

JR+ Jo Jo 



On the other hand, we apply a first discrete integration by part with respect to i G {0, . . . , n^ — 



1} to the second term which can be written as 



A jy- /rpn+l\7/2 _ /„n+l^'i'/2 

4 ^ -ijiA.EE 'H.Jl.' (''^■-^? 

neN i=0 * 

- 7 At E - 7 At E T:X\ <-i- 

neN neN 

Then, introducing DhTh a discrete approximation of the gradient of Th by 

"3~1 rpn r-pn 

DhTh{t,s) := E E^ A^+Ag' i ^{s.,s,+i)(') Mt^,t^+n{t), 
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we have 

2K\\ f /-Sn.-l , 

£l = ^ / DhTl^\t + At,s)dsip{t,s)dsdt 
I Jr+ J so 

-if Thit + At,0)ip{t,so)+Thit + At,l)ip{t,Sn,^i)dt 
Passing to the hmit /i ^ 0, we get that 

£l^^^ [ [ dsT'^/^{t,s)ds^{t,s)dsdt-i [ T{t,0)ip{t,0) +T{t,lMt,l)dt. 
' Jr+ Jo Jr+ 

Finally, we conclude that T is a weak solution of (2.1). By uniqueness of the solution to (2.1), 
it yields that the sequence (Th)h converges to the weak solution of (2.1). 

The implicit scheme (2.10) is unconditionally stable, but it requires the numerical reso- 
lution of a nonlinear system. For this purpose a Newton method is applied which increases 
considerably the computational cost and makes this method inefficient. Another strategy 
would consist in applying a semi- implicit scheme for the time discretization, but it still re- 
quires the implementation of a new linear system at each time iteration and the computational 
cost remains too important. In the following we propose a numerical scheme inspired by the 
work of F. Filbet & S. Jin [ ] to handle with this problem. 



2.6. An implicit-explicit (IMEX) scheme. In [ ], the authors proposed to handle with a 
stiff and nonlinear problem. The main point is to write the nonlinear problem in a different 
form in order to split the nonlinear operator in the sum of a dissipative linear part, which can 
be solved in an implicit way and a non dissipative and nonlinear part which will be solved 
with a time explicit solver. The main difficulty is to find an adequate decomposition of the 
operator. For instance the nonlinear diffusive operator can be written as 



dlT + ds Kn r5/2 _ ^] Q T 



and the time discretization to (2.1) becomes 



(2.15) 



At 



((^11 



■n\5/2 



z.a,r"+i(o) + 7r"+^(o) = [K\\ (r"(o))'/" - a,r"(o) 



z^a,r"+i(i) -7r"+^(i) = [k» (r"(i))°/^ -z^) a,r"(i). 



5/2 



To choose an appropriate v for the scheme (2.15), we perform an energy estimate of the 
numerical approximation. 

Proposition 2.6. Assume that the viscosity term u is such that 

(2.16) K\\ ||r"||^^ < V, Vn G N. 
Then the numerical solution satisfies the following 

(2.17) \j\T-+^ fds + ^ j\dsT-+^\'ds < \j\T-fds + ^ j\dsT-\' ds. 
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Proof. We multiply (2.15) by T^~^^ and integrate on s € (0, 1), hence we have 



^l|2 



< 



J 



(is 



A*7((To"+^)%(t:;_\)1- 



1 \2 



Using the assumption that \^n\^l'^ ^ ^ g^j^^ applying the Young's inequality, we obtain 

2 



V - Kii T' 



1 5/2 



2e 



Therefore with the choice e = i^, we have 



- C {T'+^f ds + - At r IdsT'^+^f ds < - C {T''fds + - At f\dsT''\^ds. 
2 Jo 2 Jo 2 Jo 2 Jo 

5 /2 

Hence, the scheme (2.15) is stable when /C||||T"|[oo < ^■ 

Now, we can give the fully discrete scheme, called in the sequel IMEX, as follows 



□ 



f , 1 T-n+1/2 

' rpn + l rpn T — T- ' I 



(2.18) 



1/2 



At 



As,; 



, < i < — 1, 



, Tf = To,i, < i < n, - 1, 



with the numerical flux J^^^y2 given for i e {0, . . . , — 2} by 



(2.19) 7;"+/// = 2 K|| 



TTft rpn rpn+l rpn+l 

^i+1 _|_ 2 i+1 i 

Asj+i + Asj 



Asj+i + Asj 



whereas at the boundary s = and s = 1, we apply the boundary conditions written in the 
form (2.15), 

-1, 



(2.20) 



n+1/2 



i+1/2 



+7^0^+1, ifi 

-.71+ 1 



I -7rr-\, if^ = n.-l. 



Moreover, the viscosity > is initially chosen as an upper bound of -fCy HT*^!!^^ and is then 
readjusted along iterations n € N in order to satisfy the condition (2.16): 



Algorithm to compute v 

V := 2i^|j||r°||^^ and n = 
while n < Nr , do 

— end 

compute the numerical solution T' 
if < |K||||r"+i||^^ then 

u ^ 2i^||||r"+i||^^ 

end if 
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if > 4K||||r"+i||^^ then 

u ^ K|il|r"+i||^V2 

end if 

n n + 1 
end while 

2.7. Numerical results. To compare the numerical results obtained with the different 
schemes, we take -) = 2, K\\ = 1 and the initial temperature is = 5, whereas the fi- 
nal time of the numerical simulation is equal to Tg^^ = 1. On the one hand a reference 
solution is computed using the finite volume method with an explicit scheme (2.8) on a uni- 
form grid with Ug = 450. On the other hand, we basically compare both implicit (2.10) and 
IMEX (2.18)-(2.20) schemes with different uniform grids with Ug = 50, 150. Furthermore, 
we choose the time step equal to At = 10^^, 10~^, 10~^ and 10~^ respectively. 



At 


10^^ 


10"^ 


10^4 


10"^ 


Implicit scheme (2.10) 


Ug = 50 


0.05 


0.31 


2.24 


22.06 


Us = 150 


0.60 


4.07 


27.49 


249.61 


IMEX scheme (2.18)-(2.20) 


Us = 50 


0.01 


0.09 


0.63 


5.34 


Us = 150 


0.10 


0.24 


2.24 


21.97 



Table 1. Computational time for the implicit scheme (2.10) and the IMEX 
scheme (2.18)-(2.20) in seconds at the final time of the numerical simulation 

Tend I- 



We observe from Table 1 that the IMEX scheme is much more efficient than the implicit 
scheme in terms of computational cost since the linear system corresponding to the implicit 
part does not depend on the iteration n when the viscosity i/ > is large enough. For Ug = 50, 
the computational time of the IMEX scheme is less than one fourth of the one corresponding 
to the implicit scheme whereas for Ug = 150, the implicit scheme is ten times more consuming 
than IMEX scheme. 



At 


10"^ 


10"^ 


10-* 


10"*^ 


Implicit scheme (2.10) 


Ug = 50 


0.0580 


0.0612 


0.0617 


0.0617 


Ug = 150 


0.0190 


0.0184 


0.0187 


0.0187 


IMEX scheme (2.18)-(2.20) 


Ug = 50 


0.0621 


0.0600 


0.0598 


0.0598 


Ug = 150 


0.0213 


0.0182 


0.0181 


0.0181 



Table 2. Relative errors obtained using an implicit scheme, IMEX scheme 
at time T^nd = 1- 



Concerning the accuracy and stability, Table 2 shows that the numerical solution computed 
with both implicit and IMEX schemes is stable for any time step At and the numerical errors 
are of the same order. Moreover, we get similar results when time step is smaller than 10"^. 
Of course, when we increase the number of points Ug, the numerical error decreases and 
the IMEX scheme (2.18)-(2.20) seems to be more accurate for small time steps. Finally, in 
Figure 2a, we observe that the large errors appear around the boundary, where large gradients 
of temperature occur. The Figure 2b illustrates the temperature evolution at different time 
t = 0.25, 0.50, 0.75 and 1. We note that the temperature has a fast decay at the beginning, 
then it stabilizes to a steady state when t approaches the final time T^nd = 1- Furthermore we 
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observe that the temperature develops steep gradients at the boundary modehng the coohng 
of the plasma due to the limiter effects. Indeed, on the one hand the thermal diffusion depends 
on the term T^l'^ which is large at the beginning and then becomes smaller and smaller. On 
the other hand, due to the nonlinear flux at the boundary when the temperature becomes 
small, the temperature gradient becomes larger and larger. 



- Reference solution 
N =50 




0.4 

s axis 



(a) 



1.6 




Time=1/4 






TimB=1/2 


1.4 




Time=3/4 


£ 1.2 




Time=1 



0.4 o.e 

s axis 



(b) 



Figure 2. Temperature evolution of problem (2.1). We use the IMEX scheme 
to approximate (2.1) and choose time step of At = 10~^. (a) reference solution 
and the results of IMEX scheme for rig = 50, 150 at time T^nd = 1, (b) results 
of IMEX scheme for = 150 at time t = 0.25, 0.50, 0.75 and 1. 



3. The 2D problem 

In this section, we consider the two dimensional problem where the temperature T depends 
on time t and two space variables (s,r) € = (0, 1) x (0,1) with appropriate boundary 
conditions 

(3.1) dtT - ds{K\\T^I^ dsT) - dr{K±drT) = 0, t > 0, (s,r) € Q, 

where i^y and K± are nonnegative constants with K± <C i^y. For the boundary conditions 
we impose a boundary flux in r = and assume that for r = 1 the flux of temperature is 
zero, that is, 

( drT{t,s,0) = -Q±, s e (0, 1), r = 0, t > 0, 

(3.2) 

[ drT{t, s, 1) = 0, s G (0, 1), r = 1, t > 0, 

and at the boundary s = and s = 1 we consider either periodic boundary conditions or of 
modelling describing the effects of the limiter which allows to decrease the temperature in 
the device. At s = 0, we have 

KuT^/^t,0,r)dsT{t,0,r) = jT{t,0,r), r G (1/2, 1), t > 0, 



(3.3) 

and s = 1, 
(3.4) 



T{t,0,r) = r(t,l,r), re (0,1/2), t>0, 

K\\T'^/^{t,l,r)d,T{t,l,r) = -jT{t,l,r), r G (1/2, 1), t > 0, 
r(t,0,r) = r(t,l,r), rG(0,l/2), t>0. 
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This model also satisfies an energy estimate given by 



2dt 



[ \T{t,s,r)\^dsdr = [ \dsT^^^\^dsdr - Ks_ [ \drT\^dsdr 

- 7 / {T{0,r) + T{l,r))dr + K^Q^ [ T{s,0)ds. 

Jl/2 Jo 



1/2 

To discretize the system (3.1)-(3.4), we apply a finite volume method in space coupled with 
a time splitting scheme for the time discretization. We first present the numerical scheme 
and describe precisely the discretization of the boundary conditions. Finally we compare our 
numerical results with those obtained by standard explicit and implicit time discretizations. 

3.1. Time splitting scheme. We apply a time splitting scheme in both directions. As for 
the one dimensional case, we apply an IMEX scheme to treat the nonlinear equation and find 
a condition on the viscosity > to get a uniformly stable scheme. We first consider the 
non linear problem in the s direction, 

(3.5) ^^^^ - ((/^||(T")5/2 - ^) a,r") - ud]x = 0, {s,r) G n, 

with the boundary condition (3.3), 

{(K\\ (r"(0,r))^/2 - u) a,r"(0,r) = 7r*(0,r) - iydsT*{0,r), r G (1/2,1), 
T%0,r)=T*{l,r), r G (0, 1/2), 
and then the condition (3.4), 

' (T"(l,r))^/2 _ \ a,T"(l,r) = -7^^(1,0 " z^5,T^(l,r), rG (1/2,1), 

(3.7) <^ ^ ^ 

^ T*il,r)=T*{0,r), r G (0, 1/2), 

which allows to compute a first approximation T*. Then we compute a numerical approxi- 
mation of the linear heat equation, 

(3.8) — dr{K^drT^+') = 0, (s,r)Gf], 

with non homogeneous Neumann boundary conditions 

a,r"+i(s,0) = -Q^, sG(0,l),r = 0, 

(3.9) 



a^r"+i(s,l) =0, s G (0,1), r = 1. 



For the sake of clarity we present a stability estimate on this semi-discrete scheme (discrete 
in time and continuous in space), but the proof can be easily adapted to the fully discrete 
case. 

Proposition 3.1. Assume that the viscosity term v is such that for any r G (0, 1), 

K\\ \\T''fJ^ < Vn G N. 
Then the numerical solution satisfies the following 

]- I {T^+^)^ drds+'^ I |9,T"+i|2drds < ]- j {T^f drds+'^ [ \dsT'>\^drds 

n+l 



K±Aty] / \drT''\'^drds - Q± / T''{s,0)ds 
fc^i Un Jo 
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Proof. Multiplying (3.5) by T* and integrating in 17, we obtain 

' d.T'' d,T* dr ds 



ly (^(^T*f -{T'^f^ drds < -At J (K||(T" 

-At [ u\dsT*\^drds 
Jn 

-7 At / |r*(0,r)|2 + |r^(l,r)|2(ir. 

Jl/2 

Then, applying the Young inequality and taking v such that for all r G (0, 1), 

< K\\ |r"(s,r)|^/^ < u, Vn G N, 



we have 



(3.10) I [ {T*f drds+ ^ / \dsT*\^drds < ]- f {T^'f drds+'^ [ IdsT'^fdrds. 
^ Jn ^ Jn ^ Jn ^ Jn 

Similarly, we multiply (3.8) by T"^-*^ and integrate with respect to (s,r) G we get 



' rpr, 



-n2 



(3.11) 



{T*f^drds < -AtK± j [dr-T'^+^f drds 

+ AtQ±K^ [ r"+i(s,0)ds. 
Jo 



Furthermore, we derive (3.8) with respect to s and get 

Then we multiply this latter equality by vdsT^^^ and integrate over (s, r) G $7, 



{dsT 



n+l 



drds < -2AtuKi 



\drsT 



n+l\2 



dr ds 



+ uAt [ds{drT^^')dsT^+']lzl. 



Hence using that ds (a,r"+i(s,r)) = 0, r G {0, 1}, it yields 



drds < 0. 



(3.12) / u\{dsT^+'f -idsT*y 

Jn L 

Then, gathering (3.11) and (3.12), we get 

1 f ^^^^ + ^ [ [|a,r"+i|2 + |9,r"+^i2] drds - aik^q^ [\''+\s,o)ds 

^ Jn ^ Jn Jo 

<- I {T^f drds+'^ I \dsT''\^drds. 

Finally, the latter inequality together with (3.10), it gives 

1 1 ^rj^n+i^ drds + ^ / [|a,r"+i|2 + E:^|a,r"+i|2] drds - AtK^Q^ [\''+\s,0)ds 
^ Jn ^ Jn Jo 

< i / (n^ drds+'-^ j |a.r«p 

By induction and summing over A; = 0, . . . , n, we get the result 

^ / (r"+i)'drds + ^ / \dsT^+^\^drds < i / {T^f drds+'^ f \dsT^\^drds 

n+l r |. /"l - 

-i^xAtV / la^r^pdrds - / T''{s,Q)ds . 

fc^i L^c Jo 



dr ds. 



NUMERICAL STUDY OF A NONLINEAR HEAT EQUATION FOR PLASMAS 



17 



□ 

3.2. A finite volume approximation. For the space discretization, we consider a set of 
points (sj_i/2)o<i<ns a set of points of the interval (0, 1) with s_i/2 = •Sns-i/2 = 1 a-^^d 
Us + l represents the number of discrete points in the direction s and ii"j-i/2)o<j<nr a set of 
points of the interval (0, 1) with r_i/2 = 0, r„^_i/2 = 1 and + 1 represents the number of 
discrete points in the direction r. For < i < — 1, < j < — 1, we define the control cell 
Ci,j by Cij = (si_i/2,Sj+i/2) X (rj_i/2,rj+i/2)- We also denote by {ri,Si) the center of Cjj 
and by Asi the space step Asi = Si+1/2 — Si-i/2 and Arj the space step Arj = rj+1/2 — rj_i/2 
where we assume that there exists ^ G (0, 1) such that 

(3.13) < Asi, Arj < h, V(ii, j) G {0, . . . , - 1} x {0, . . . , - 1}, 

with h = maxjj{Asj, Arj}. 

We shall construct a set of approximations Tij{t) of the average of the solution to (1-2)- 
(1.3) on the control volume Cij and set 



Hence, the finite volume discretization to (3.5) can be written as 

T* -T" jr"+l/2 _^n+l/2 

= -^/^^ V(^,J)G{0,...,n.-l}x{0,...,n.-l}, 

At Asi 

where the flux J^i-\-i/2j corresponds to the one dimensional flux given by (2.19) and periodic 
boundary conditions are applied for rj G (0, 1/2) and conditions (2.20) for rj G (1/2, 1). 
Then, the finite volume discretization to (3.5) can be written as 



rjin + l rpi, _ 

At Avj 



G {0,...,n, -1} X {0,...,n^-l} 



where Gij+1/2 is given by 

rj-m + l rpn+l 

(3.14) ^^,,+1/2 = 2 ^i;' , j = 0,...,nr-2. 

' Arjj^i + Arj 

Moreover, at the boundary r = and r = 1, we apply the boundary conditions, 

f -K^Q^, ifj = -l, 

(3.15) C?ij+i/2 = < 

[ 0, if j = nr- 1. 

3.3. Numerical results. In this section we compare the different numerical results related 
to the 2D problem (3.1)-(3.4) obtained using a time splitting scheme with an explicit, implicit 
and IMEX treatment of each step. As before, we first compute a reference solutions obtained 
from an explicit scheme with a small time step satisfying a CFL condition At ~ /i^. In 
the following numerical simulations, we choose the different physical parameters as K^^ = 1, 
K± = 10~^, 7 = 2, Q± = 10. Moreover, the initial temperature is given by 

(3.16) r0(s,r)=3, 

and the final time of the simulation is T^^nd = 2. 

To compute the reference solution, we have chosen rig = 300 and = 300, whereas the 
numerical results using implicit and IMEX schemes are obtained with Us = 100 and rir = 100 
with several time steps At = 10~^, 10~^, 10~^, and 10~^. First, concerning the computational 
time we observe in Table 3, that the IMEX scheme is much faster than the implicit scheme. 
Furthermore, the numerical error presented in Table 4 for both scheme is of the same order of 
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magnitude and thus the IMEX scheme is clearly much more efficient than the fully implicit 
scheme. 



At 


10-^ 


10"^ 


10"^ 


10-^ 


Implicit scheme 


4.02 


25.64 


172.95 


1327.50 


IMEX scheme 


1.62 


4.42 


36.24 


403.63 



Table 3. Computational time for the 2D problem (3.1)-(3.4) using implicit 
and IMEX schemes at time Tf,nd = 2. 



At 


10-^ 


10"^ 


10~^ 


10-^ 


Implicit scheme 


0.2245 


0.0236 


0.0020 


2.1985e-04 


IMEX scheme 


0.2093 


0.0213 


0.0018 


2.4385e-04 



Table 4. The relative errors for the implicit and IMEX schemes compared 
with a reference solution for the 2D problem (3.1)-(3.4) at time Tend = 2. 

Now we want to investigate the effect of the splitting scheme on the numerical error and the 
computational cost. Therefore, we also propose a comparison between the different schemes. 
We first compare the computational time applying the IMEX scheme with and without the 
splitting method with a time step At = lO'^, (n^.n^) = (50,50), (100,100), (300,300) and 
(500, 500) respectively. On the one hand, we observe in Table 5 that the splitting method is 
much faster than the non-splitting method when the number of discrete points increases. 



71/ g X Tt'p 


50 X 50 


100 X 100 


300 X 300 


500 X 500 


IMEX Non-splitting scheme 


11 


60 


505 


2112 


IMEX splitting scheme 


16 


36 


219 


601 



Table 5. Computational time of IMEX with and without splitting scheme 
at time T^nd = 2. 



On the other hand, we compare the numerical errors corresponding to the two strategies 
with {ns,nr) = (100, 100), At = 10~^ in Table 6, in particular the fully implicit scheme with 
and without splitting and the IMEX scheme with and without splitting. We observe that the 
method without splitting is always more accurate than the one with the splitting method. 



Scheme 


Splitting implicit 


Splitting IMEX 


Implicit 


IMEX 


Numerical error 


2. X 10~^ 


2. X 10"^ 


5. X lO-"* 


5. X lO""* 



Table 6. Relative errors for different numerical schemes compared with a 
reference solution for (n<j,nr) = (100, 100), At = 10~'^ at time T^nd = 2. 



In Figure 3, we present the evolution of the approximation of the temperature (3.1)-(3.4) 
in computational domain Q, which is divided into two regions : the transition layer and the 
scrape-off layer (SOL) as illustrated in Figure 1. We first initialize the temperature to a 
constant and then observe immediately that temperature decreases rapidly in the scrape-off 
layer and becomes singular around the limiter (which corresponds to the boundary s = and 
1 with r > 1/2). On the other hand, in the transition layer, the temperature converges to 
a steady state which is homogeneous in s € (0, 1). The different numerical schemes give the 
same qualitative behavior of the solution. 
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{e)t = l (f) i = 2 

Figure 3. Temperature evolution of problem (3.1). 



2.245 



1.194 



0.5092 



0.2 0.4 0.6 0.8 



Max: 4.518 
Min: 0.01818 



In Figure 4, we plot the temperature evolution at the section r = 0.25, r = 0.75 and 
s = 10~^ and s = 0.5 respectively. According to Kocan et al. [13, 14], the parallel thermal 
diffusivity is much larger than the perpendicular one, i.e. K^^ ^ K±. Therefore, the tem- 
perature becomes constant along the magnetic field lines, that is for s € (0, 1). We observe 
in Figures 4 that the temperature is constant at all time whereas steep gradients develop 
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at the boundary layer s = and s = 1) in the SOL region. In the perpendicular direction 
r, the situation is different. We also observe that at time t = 2 the temperature decreases 
linearly with respect to r in the transition layer (0 < r < 0.5), according to the heat flux Q± 
at edge r = 0, and then decreases exponentially in the scrape-off layer (0.5 < r < 1). These 
numerical results correspond to the retarding field analyzer (RFA) [12, 13, 14]. 
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Figure 4. Temperature evolution at section r = 1/4, r = 3/4, s = 10 ^ and 
s = 1/2 at time t = 0.1, 0.5, 1 and 2 respectively. 



Finally, we present the evolution of the energy dissipation with respect to time: 



ld_ 
2di 



I \T{t,s,r)\'^dsdr = £i + £2 + £3, 
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with 

£i := - f (a'ii (r)^/2 \dsT\^ + K± \drT\^^ drds, 
< £2 ■■= -7 / T{t,l,rf + T{t,0,r)^ dr, 

Jl/2 

£3 := +Q±K^ [ T{s,0)ds. 

V JO 

The Figure 5 states the terms fi, £2, £3 as function of t. We plot these terms obtained from 
imphcit and IMEX schemes. Note that these two figures are almost the same. In fact, at 
the beginning of simulation, there is a fast decay of the temperature, thus the quantity —£i 
representing the total energy exchange ratio in the domain is increasing for t < 0.1. Then, 
it converges to an equilibrium state for larger time. On the other hand, the quantity —£2 
decreases with respect to time, it is due to the anisotropy between K^^ and Kj_. Indeed, the 
energy is transferred to the limiters in the scrape-off layer region whereas in the perpendicular 
direction r, the thermal diffusivity is small. Finally, as we have seen in Figure 4 on the edge 
of of the core, the temperature does not vary significantly, thus the quantity £3 increases 
slightly with respect to time. 




Time Time 



(a) Implicit scheme (b) IMEX scheme 

Figure 5 . Evolution of the energy dissipation with respect to time for prob- 
lem (3.1), with At = 0.001. 



4. The coupling problem 

In this section, we consider the full 2D model (1.2) composed of two different particle 
species, i.e. ions and electrons. We denote by Tj {resp. T^) the temperature of ions [resp. 
electrons) which depends on time t and two space variables (s, r) G O. The two equations are 
coupled by a non-zero source term which balances the temperature between the two particle 
species, 

r dtTi - ds{K\\,iT^'^dsTi) - driK^^idrTi) = +P{Ti - Te), for {s,r) G Q, 

(4.1) 

[ dtTe - ds{K\\^eT!'^dsTe) - driK^^ATe) = -m " Te), foT {s, r) € n, 
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where K±^i <C K± e <C -fCn g and /3 is a negative constant. These two equations are 

completed with the same type of boundary conditions as in (3.2)-(3.4). 

4.1. Time splitting scheme. Now we discretize the full system (4.1) using a splitting 
scheme in three steps. We assume that an approximation of the solution {Te,Ti) at time 
is known and denote it by {T",T^). Therefore, we first approximate the source part 
coupling the two temperatures and Tj using an implicit scheme, which yields 



(4.2) 



T* = - (l 1 T" + - ( 1 H 1 

T* = -(1 H )r" + -fl )T" 

2^^ 1 -2/3At^^ ^ 2^ l-2/3At^ * 



It is clear that (4.2) guarantees the positivity of the temperature. Then we apply the same 
time splitting steps as before in direction s and in direction r as follows. On the one hand 
we compute T** for a G {i, e} by solving (3.5)-(3.7). On the other hand we apply the last 
step (3.8)-(3.9) in the direction r. 

Furthermore, for the scheme (3.5)-(3.7), (3.8)-(3.9) and (4.2), we also prove an energy 
estimate 

Proposition 4.1. Consider that the initial datum Tq is nonnegative and Tq G L°°(0, 1). 
Assume that for a G {i, e}, the viscosity term v is such that for any r G (0, 1), 

max i^ll^JlT^ 11^2 < y^^^^ 

Then the numerical solution, given by (4--2), satisfies the following 

\ E / {\TT'? ^ ^tv%T-'^\^\drds 



< ^ E / {\Tl? ^ ^tv\d,Tl\^\drds 

n+1 



- At ^ E^^." / \QrTlf'^rds 

n+1 »i 

+ At E^^'"^^'" / Ta(S'0)ds. 



i£{i,e} ; 

Proof. We first observe that the energy estimate of the two last steps in the direction s and 
r are the same as the one proved in Proposition 3.1, hence we have 

i f [iT^ip + uAt (la.rrT + Kx,„|9,rri|2)] drds 
^ Jn 

< \ ( [\Tt:? + Atz^ia^rrHdrds + Ati^^,«Qx,a / T^+\s,^)ds. 
^ Jn Jo 

Therefore, to achieve the proof on the energy estimate, we only observe that (4.2) can be 
written as follows 



(4.3) 



-T^ = +At/3 {Tt - Tt) 
- T2 = - At/3 {T* - 



e . 
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Multiplying the first equation (4.3) by T* and the second by T* and integrating on (r, s) G 
it yields 

I- [ \T^\^ + \T^\^drds <]- ! \Tr\^ + iTJ^pdrds 
2 ./o 2 



- E 

2 ^ 



Moreover, differentiating (4.3) with respect to s and multiplying the first equation by udgT* 
and the second one by udgT*, we get 

^ / %T:\^ + \dsTt\^dTds <^ [ \d,T.r\' + \dsT^\^drds. 
^ Jo. ^ Jn 

Finally, we have 

/ [|T„"+^P +At (z.|9,r^+i|2 + Kx,a|5.r;:+i|2)] drds 
Jn 

< \Y. j +^t[^\dsK? + K^^^\drT^\^)]drds 

o£{i,e} ^ 

- At ^ K±,ag±,a rs^\s,o)ds. 

Summing over k = 0, . . . ,n, we complete the proof. □ 

Finally space discretization is performed using the finite volume scheme presented in Sec- 
tion 3.2. 

4.2. Numerical results. In this section, we compare the numerical results obtained from 
the implicit scheme and the IMEX scheme for (4.1). We choose j = 2 x 0.01, i^'y g = 1, 
= 0.01, K^^e = 0.01, 7j = 0, 7e = 2.5, = = 10 and (3 = -0.02. The' initial 
temperature is such that 

T^{s,r) = 3, and T°{s,r) = 3, (s,r) € fl. 

The final time of the simulation is T^^d = 1 and the mesh size is chosen as Ug = 100, 
rir = 100. 

We plot the electron and ion temperature and compare their ratio at different time. The 
aim is to compare the different behaviors between electron and ion temperatures at the edges 
and in the scrape-off layer of a Tokamak ['. ]. 

On the one hand, we propose in Figure 6, the temperature evolution. On the left hand 
side, we present the electron temperature, whereas on the right hand side we give the ion 
temperature. We first notice that the electron parallel thermal diffusivity is about 100 times 
larger than the one for ions [2, 10], and the electron energy exchange ratio at the edge 

r € (0.5, 1) depends on 0{Te ), thus the temperature has a fast decay when it is small 
in the scrape-off layer. However, the boundary conditions for ions in the scrape-off layer is 
given by the homogeneous Neumann condition dgTi = 0, which means that there is no energy 
exchange at the limiters. Thus the ion temperature does not vary significantly at scrape-off 
layer. 



On the other hand, the ratio between electron temperature and ion temperature is pre- 
sented in Figure 7. The Figure 7 illustrates that in the transition layer, the ion and electron 
temperatures are almost identical. However, in the scrape-off layer, at the final time T(,nd = 1 
the ratio r becomes large around the limiters due to the boundary condition dsTg cc Tg . 
The evolution of the ratio r in the radial direction is given in Figures 7. We observe that in 
the transition layer the ratio r is almost equal to 1, whereas in the scrape-off layer this ratio 
becomes large. For example, at time t = 1 the ratio r = 6 for s = 1/2 while it is r = 45 for 
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(b) Tj at t = 0.25 




(d) Ti at t = 0.5 
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(d) Te at t = 1 



(e) Ti at t = 1 



Figure 6. Temperature evolution of problem (4.1). 



s = 10~^. These behaviors correspond to the experiment results in Kocan et al. [1.3, 14]. At 
last we vary the parameter (3 to study the equilibrium source term in Figure 8 and observe 
that when the parameter |/3| is large, the ratio r decreases. 
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Figure 7. Ratio r = TijT^ at section r = 1/4, r = 3/4, s = 10"^ and s = 1/2 
at time t = 0.1, 0.25 and 1 respectively. 
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Figure 8. Ratio r = Tj/Tg at section s = 10 ^ and s = 1/2 for different 
parameters /3 = — 2 x 10"^, —2 x 10^^ and —2 at time t = 1. 
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5. Conclusion 



We have presented various numerical approximations for a nonlinear temperatTirc balance 
equation describing the heat evolution of a magnetically confined plasma in the edge region 
of a tokamak. Numerical comparisons show that an IMEX scheme based on a "smart" 
decomposition of the nonlinear diffusive operator coupled with a splitting strategy gives an 
efficient numerical scheme in terms of accuracy, stability and reasonable computational cost. 
The next step would consists to couple the present model with the transport equations for 
the plasma density and momentum. 
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